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Abstract. We examine the phenomenon of dynamical heterogeneity in computer 
simulations of an equilibrium, glass-forming liquid. We describe several approaches 
to quantify the spatial correlation of single-particle motion, and show that 
spatial correlations of particle displacements become increasingly long-range as the 
temperature decreases toward the mode coupling critical temperature. 
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1. Introduction 

Liquids cooled toward their glass transition exhibit remarkable dynamical behavior 
The initial slowing down of many liquids at temperatures well above their glass transition 
temperature T g can be described to great extent by the mode coupling theory |2], ]3|, §], 
which predicts diverging relaxation times at a dynamical critical temperature T c despite 
the absence of a diverging or even growing static correlation length ||. At the same 
time, experiments in the temperature range T g < T < T c , and simulations in the range 
T > T c , have shown that it is possible to select subsets of particles in the liquid (or 
monomers in the case of polymers) that move differently from the bulk on time scales 
less than the structural relaxation time H 0, H 0- The question then arises as to 
whether glass-forming liquids exhibit spatially heterogeneous dynamics, and if so, is 
there a growing length scale associated with this dynamical heterogeneity? 

We have employed three complementary computational approaches to address this 
question: 

(i) Subset approach. In this approach ]T0L [TTJ, |T^] , we monitor the displacement of each 
particle in a time window At, rank the displacements from largest to smallest, and 
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examine spatial correlations within subsets of particles exhibiting either extremely 
large or extremely small displacements. In this way, we have demonstrated that 
particles of similar mobility form clusters that grow with decreasing temperature 
T. Moreover, particles exhibiting the largest displacements within a time window 
in the late-/?/early-a relaxation regime ("mobile" particles) form low-dimensional 
clusters that grow and percolate at the mode coupling transition temperature T c . 



(ii) Displacement- displacement correlation function approach. In this approach |13|, |14|, 
[15|| , we monitor the displacement of each particle in a time window At as above, 
and calculate a bulk, equilibrium correlation function that measures the correlations 
between the local fluctuations in the particle displacements. Because this approach 
does not rely on defining subsets via arbitrary thresholds or cutoffs, it may prove 
more amenable to future analytical treatment. 

(iii) First passage time approach. In this approach [|l(J, we monitor the first passage 
time r for each particle to move a distance e, and then for a given e, examine spatial 
correlations between the different r values either by grouping particles into subsets 
of similar r, or by calculating a r-r correlation function using all the particles. This 
approach shows that particles with similar first passage times to travel any e are 
spatially correlated. It also allows us to identify particular values of e for which the 
first passage times are most correlated, and thereby extract detailed information 
on, e.g., the cage size. 

In this paper, we present results from the first two approaches outlined above, and 
show that as an equilibrium Lennard- Jones liquid approaches T c , correlations between 
local fluctuations in particle displacements develop, grow, become increasingly long- 
ranged, and appear to diverge at T c . Similar results for a simulated polymer blend are 



presented elsewhere [fL7 |. 



2. Theory 

Consider a liquid confined to a volume V, consisting of N identical particles, each with 
no internal degrees of freedom. Let the position of each particle i be denoted T*j. The 
motion of a particle i may be described by its displacement /^(t, At) = |rj(t+At) — Tj(t)| 
over some interval of time At, starting from a reference time t. The question of 
whether particle motions in a liquid are correlated can be addressed by modifying 
the usual definition of the density-density correlation function G(r) (which measures 



spatial correlations of fluctuations in local density away from the average value [[18], [L9| ) 
so that the contribution of each particle i to G(r) is weighted by /ij. Thus, we define a 
"displacement-displacement" correlation function |L3|, [14|, [15], f20| , 

G u (r, At) = J dr'([u(r' + r,t,At) - (u)\ [u(r',t,At) - (u)\), (1) 

where 

N 

u{r, t, At) = £ fjn(t, At) 5{r - r»(t)). (2) 
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and where (• • •) denotes an ensemble average. G u (r,At) measures the correlations 
of fluctuations of local displacements away from their average value. Since we are 
considering an equilibrium liquid, and because G u is defined as an ensemble average, 
G u does not depend on the choice of the reference time t. For the same reasons, 
(u) = (u(r,t, At)) only depends on At for a homogeneous, equilibrium liquid. 

G u (r, At) can be separated into self and distinct parts, so that we can identify a 
spatial correlation function g u (r, At) analogous to the static pair correlation function 
g{r) conventionally used to characterize the structure of a liquid jl8| : 



G u (r,At) = (N)^S(r) + (u)(U) \g u (r,At) - 1], (3) 

where 

Y INN \ 

9u(r, At) = ^7y(E E M*, At ) Pitt, At ) S ( r + ~ »"*(*)) )• (4) 



=1 3=1 
3^i 



Here we have defined the "total displacement" U as U(t,At) = f dru(r,t,At), and its 
ensemble average (U) = (U(t,At)). We have also defined p, = (jr J2iLi Pi{t, At)), and 
A* 2 = {jt Sili lA {t, At)). In a constant-iV ensemble, (u) and (U) are readily evaluated 
as (u) = pp. and (U) = N ft. Note that g«(r, At) is different from the distinct part of the 
time-dependent van Hove correlation function g(r,t), because it correlates information 
between pairs of particles using information about the position of each particle at two 
different times. Thus g u (r, At) is in some sense a four-point correlation function, whereas 
g{r,t) is a two-point correlation function. 

A "structure factor" for the particle displacements can be defined as 

/ 1 N N \ 
S u (q, At) = ( —= £ £ tk(t, At)tij(t, At) exp [-iq ■ (r^t) - rj(t))} ).(5) 
\Nfi 2 i=lj=1 I 



The fluctuations of U are related to the volume integral of G u (r, At): 

([U-(U)] 2 ) = J drG u (r,At) = (U)(u)kTK u , (6) 

where k is Boltzmann's constant. Recall that the volume integral of the density- 
density correlation function G(r) is proportional to the isothermal compressibility k, 
which diverges at a conventional critical point because G(r) becomes long-ranged. 
Likewise, the fluctuations in U will provide information regarding the range of G u (r, At). 
Thus, in analogy with conventional critical phenomena, we have defined an isothermal 
"susceptibility" k u . 

For a fixed choice of At, note that if the displacement were always the same for 
every particle, then g u (r,At) and g{r) would be identical for all r (and S u (q, At) and 
S(q) would be identical for all q). Hence, it is deviations of g u (r,At) (S u (q,At)) from 
g(r) {S(q)) that will inform us of correlations of fluctuations of local displacements 
away from the average value, that are in excess of those that would be expected based 
on a knowledge of g{r) or S(q) alone. 
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3. Simulation Details 

We have measured the spatial correlations in particle displacements using data obtained 
PH |TT| , O] from a molecular dynamics simulation of a model Lennard- Jones glass- 
former. The system is a three-dimensional binary mixture (80:20) of 8000 particles 
interacting via Lennard- Jones interaction parameters pl| . We analyze data from seven 



(p, P, T) state points that lie on a line in the P, T plane, approaching the mode- 
coupling dynamical critical temperature T c « 0.435 at a pressure P c ~ 3.03 and 



density p c ~ 1.2 |22| (all values are quoted in reduced units f21|l ). The highest and 
lowest temperature state points simulated are (p = 1.09, P = 0.50, T = 0.550) and 
(p = 1.19, P = 2.68, T = 0.451), respectively. Following equilibration at each state 
point, the particle displacements are monitored in the NVE ensemble for up to 1.2 x 10 4 
Lennard- Jones time units (25.4 ns in argon units) for the coldest T. Complete simulation 
details may be found in |I2| . 

For all seven state points, a "plateau" exists in both the mean square displacement 
and the intermediate scattering function F s (q,t) as a function of t fl2"fl . The plateau 
separates an early time ballistic regime from a late time diffusive regime, and indicates 
"caging" of the particles typical of low T, high density liquids. The same model liquid, 
simulated along a different path toward the same mode coupling critical point, was found 
to be well described by the ideal mode coupling theory |2"2|| . For the present (larger) 
simulation, the a- relaxation time r a that describes the long-time decay of F s (q, t) at the 
value of q corresponding to the first peak in the static structure factor S(q) increases 
by 2.4 orders of magnitude, and follows a power law r a ~ (T — T c )~ 7 , with T c ~ 0.435 
and 7 ~ 2.8. The diffusion coefficient D is found to behave as D ~ (T — 0.435) 2 ' 13 . 
The simulated liquid states analyzed here therefore exhibit the complex bulk relaxation 
behavior characteristic of a supercooled liquid approaching its glass transition. 

4. Results 

Throughout the following subsection, all quantities are calculated using all 8000 particles 
in the system. Results do not change substantially if the minority particles are not 
included in the analysis. In the cluster analysis in Subsection only the majority 
particles are included in the analysis. However, our results do not change quantitatively 
when the minority particles are included. Recall that in the mode coupling analysis of 
this model liquid by Kob and Andersen, the same mode coupling transition temperature 
was found for the majority particles as for the minority particles |22 |. 



4-1. Displacement- displacement correlation function 

S u (q,At) is shown for different T in Fig. [I], calculated for a value of At chosen on 
the order of r a for each T. For intermediate and large q, S u (q,At) coincides with the 
static structure factor S(q) = YliLi J2f=i exp [—iq ■ (r, — (see inset). However, 
for q — ► a peak in S u (q, At) develops and grows with decreasing T, suggesting the 
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Figure 1. Displacement "structure factor" S u (q, At) for different T. The q = values 
are obtained from S u (q — 0, At) = (p 2 / /i 2 ) pkTn u . INSET: Static structure factor 
S(q) for four different T. 



presence of increasing long range correlations in G u (r, At). No such growing peak at 
q = appears in the static structure factor S(q) (cf. inset), indicating the absence of 
long range correlations in the particle positions. Indeed, a comparison of g u (r, At) and 
g(r) shows that at every T simulated, the particle displacements are more correlated 
than their positions. At the lowest T, this excess correlation persists for values of r up 
to 6 interparticle distances [ [TBI . 

In Fig.|l] At is chosen on the timescale of the a-relaxation time. The significance 
of the time interval At can be seen by examining the "susceptibility" K u (At). Here 
we have calculated n u (At) directly from the fluctuations in the total displacement U . 
Fig. H demonstrates four important points: (i) As At — > 0, n u — > 0. Note that in the 
limit At — > 0, fii(t,At)/At is equal to the magnitude of the instantaneous velocity. 
Thus correlations cannot be observed by looking at a "snapshot" of the system, that 
is, by measuring correlations in the instantaneous velocity f23fl . (ii) Initially, as At 
increases, k u increases. Thus the spatial correlation of particle displacements develops 
over time, (iii) There is a time window At* during which k u is maximum. Both the 
maximum value of k u and corresponding time window At* increase with decreasing T. 
(iv) As At — > oo, n u decreases. Thus in the diffusive regime, where the particles act like 
Brownian particles, there are no correlations in the particle motions. 

Fig. ^ shows that the time At* at which the fluctuations in the particle 
displacements are a maximum follows a power law with T: an excellent fit to the form 
At* ~ (T - T c )~ 7 is obtained when T c = 0.435 ± 0.005, and yields 7 = 2.3 ± 0.3, where 
the highest value of 7 is obtained with the lower bound on T c . Our estimated value 
for 7 differs from that describing the divergence at T c of the structural relaxation time 
r a (7 ~ 2.8), but (within our numerical uncertainty) cannot be distinguished from the 
exponent governing the apparent vanishing of the diffusion coefficient D (7 ~ 2.1 ±0.1). 
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Figure 2. K u (At) as a function of At for different temperatures. 



500.0 



400.0 



300.0 



200.0 



100.0 



0.0 



10° 



10' 



10 



10" 




10 



10 X 10 



At =0.03(T-0.435)~ 



10 




0.44 0.46 0.48 0.50 _ 0.52 0.54 0.56 0.58 



T 



Figure 3. Time when the fluctuations in the particle displacements are largest, plotted 
vs. temperature. The solid curve is a fit to At* ~ (T- 0.435) ~ 23 . fNSET: At* vs. T a , 
plotted log- log. The solid line has slope 0.8f , and represents the best fit curve through 
the data. The dashed line has slope 1.0, and is not nearly as good a fit. That At* 
does not scale linearly with r Q is further supported by the fact that an exponent less 
than 2.6 describing the apparent divergence of r a with T — T c can be excluded. 



In Fig f| we show the T-dependence of n u (At*). We find that K u (At*) grows 
monotonically with decreasing T, indicating that the range of the correlation measured 
by g u (r, At*) is growing with decreasing T. The data can be fitted extremely well with 
a power law K u (At*) ~ (T — T c )~ 7 with T c = 0.435, and yields 7 = 0.84. (As shown in 
the figure, a larger (closer to Curie- Weiss-type) exponent can be obtained by fitting the 
data with a power law that diverges at a slightly lower value of T c .) Thus the previously 
determined mode coupling critical temperature coincides, within our numerical error, 
with the temperature where G u (r, At) becomes long ranged. 
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Figure 4. "Susceptibility" n u (At*) versus temperature. The solid and dashed curves 
are power law fits to the data, as indicated in the figure. 




Figure 5. Particles exhibiting the 5% largest (light) and 5% smallest (dark) 
displacements in At* starting from an arbitrary time origin, at T — 0.451. The 
clustering of particles of similar mobility is evident in the figure. 



The displacement-displacement correlation function measures the tendency for 
particles of similar mobility to be spatially correlated (cf. Fig. [5]). To extract a 
correlation length from the data, we have attempted to fit S u (q, At*) using an Orstein- 
Zernike form, S u (q, At*) oc 1/(1 + £ 2 q 2 ). This form assumes that G u (r,At*) is 
asymptotically proportional to e~ r ^/r, where £ is the correlation length (T^|. We find 
that this form fits well at the highest T but fails to fit the data on approaching T c . The 
data can instead be fitted at all T using the more general form S u (q, At*) oc 1/ (l + i^q) 11 )- 
The best fit gives r] = 1.99 and 3.54 for our highest and lowest T, respectively. However, 
because rj depends on T, £ does not show a definite trend with temperature, giving values 
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of roughly three particle diameters. Larger simulations may be required to accurately 
determine the correct functional form for S u at small q. Nevertheless, the data show 
unambiguously that as T — > T c , spatial correlations between the displacements of 
particles arise and grow and become long-ranged at T c . To obtain further evidence 
for a growing dynamical length scale, we turn now to the subset method described in 
the Introduction 



4-2. Cluster analysis of mobile particles 

Examination of the distribution of particle displacements (the self part of the usual time- 



dependent van Hove correlation function [18]) immediately shows the problems inherent 



in defining subsets of particles of extremely low mobility. Over the time range coincident 
with the plateau in the mean square displacement, this function exhibits a narrow peak 
containing most of the particles, and a long tail to large distances containing a small 
percentage of the particles. The non-Gaussian parameter 0:2 (i) has been used to identify 
a time when this tail is most pronounced [ID]. At that time, approximately 95% of the 



particles are contained in the peak, and about 5% are contained in the tail, regardless of 
temperature. Because of this tail, the most mobile particles in the liquid "distinguish" 
themselves from the bulk, while it is more difficult to define a threshold with which to 
identify particles exhibiting extremely small displacements ("immobile" particles). 

By examining the maximum displacement achieved by a particle in a time window 
At, a somewhat better definition of immobility can be achieved. It was shown in Ref. [|l^ 



that particles exhibiting the smallest displacement defined in this way are spatially 
correlated, and form relatively compact clusters. Using this definition of maximum 
displacement, which for the most mobile particles gives the same results as the simple 
definition of displacement used above, we find that highly mobile particles also cluster. 
These clusters are very ramified, and are composed of smaller "strings" of particles that 
follow one another (cf. Fig. 0). 

In this subsection, we are interested in the growth of clusters of highly mobile 
particles with decreasing temperature. To this end, we have constructed clusters of 
nearest neighbor particles whose maximum displacement in a time window At ~ At* 
falls within the top 3%, 5% and 7% of all displacements [T2). The distribution of 
clusters of size n constructed using the top 5% of the particles is shown in the inset 
of Fig. ^ for four different T. Also shown in Fig. [?] is the T-dependence of the mean 
cluster size S = J2 n2 P( n )/ J2 n P( n ) f° r eac h subset [p6fl . We find that for each subset 
S ~ (T - T p )-t with T p = 0.440,0.431, and 0.428, all close to T c . (Note that as the 
number of particles contained in the subset increases, T p decreases.) For the subset 
containing 5%, the data fall on a straight line when plotted log-log against T — T c . 
Thus, we conclude that T c appears to coincide with a percolation transition of the most 
mobile particles in the liquid. 

The presence of clusters of mobile particles whose size grows with decreasing 
temperature contributes to the growing range of the displacement-displacement 



Quantifying spatially heterogeneous dynamics in computer simulations of glass-forming liquids^ 




Figure 6. A large cluster of mobile particles found at T = 0.451, at two different 
times. Pink denotes the particles at t = and blue denotes the particles at time 
At later, where At is slightly less than At*. The small black lines connect the same 
particle at the two times. Thus between t = and At, each particle shown has moved 
almost one interparticle distance. This cluster can be decomposed into a number of 
elementary "cooperatively rearranging strings" , defined as minimal groups of particles 
that follow one another [fill Ga. 




Figure 7. Mean cluster size S plotted versus T, for subsets containing 3% (circles), 
5% (squares), and 7% (diamonds) of the most mobile particles at T — 0.451. The lines 
are power law fits 5" ~ (T-T p )-f. Best fit parameters are T p = 0.440, 0.431 and 0.428, 
respectively, and 7 = 0.397,0.687, and 0.741, respectively. INSET: Distribution of the 
size n of clusters of mobile particles for four different T. 
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correlation function. Thus we can use the mean cluster size to give a rough estimate 
of the length scale on which particle motions are correlated. In our coldest simulation, 
the mean cluster size is approximately 10 particles. Since the clusters have a fractal 
dimension of approximately 1.75 ||12|| , that implies a mean radius of gyration larger 
than 3 particle diameters. The largest cluster observed at T = 0.451 contains over 100 
particles, and thus has a radius of gyration that exceeds the size of our simulation box 
at that T. 

It is important to note that while the mean size of clusters of mobile particles 
appears to diverge at T c , the smaller, cooperatively rearranging "strings" that make 
up these clusters also grow with decreasing T, but are still finite at T c jLl]]. In the 
temperature regime studied here, the distribution of string lengths is exponential, and 
the mean string length appears to diverge at a temperature substantially lower than 
T c . Thus while the mode coupling transition may coincide with a percolation transition 
of clusters of highly mobile particles, our data is compatible with the idea of Adam 
and Gibbs that the ideal glass transition at T Q < T c is associated with the growth and 
possible divergence of the minimum size of a cooperatively rearranging region ( "strings" 
in the present simulation). This idea will be further explored elsewhere |[25|| . 

5. Discussion 

In this paper, we have defined a bulk correlation function that quantifies the spatial 
correlation of single-particle displacements in a liquid, and we have examined spatial 
correlations within subsets of highly mobile particles. While the first approach has 
the advantage of allowing a direct calculation of spatial correlations in local particle 
motions without having to define arbitrary subsets, the second approach has the 
advantage of permitting a geometrical analysis of clusters of highly mobile and highly 
immobile particles. Using both approaches, we have shown in computer simulations of 
an equilibrium Lennard- Jones liquid that the displacements of particles are spatially 
correlated |27]] over a range and time scale that both grow with decreasing T as 



the mode coupling transition is approached. While mode coupling theory currently 
makes no predictions concerning a growing dynamical correlation length, calculation 
of the vector analog of the displacement-displacement correlation function should be 
tractable within the mode-coupling framework. Finally, the displacement-displacement 
correlation function has allowed us to identify a fluctuating dynamical variable U whose 
fluctuations become longer ranged and appear to diverge at T c . In this regard, U is 
behaving much like a static order parameter on approaching a second-order static critical 
point, suggesting the possibility that we can obtain insights into the nature of glass- 
forming liquids using an extension to dynamically-defined quantities of the framework 
of ordinary critical phenomena. 



6. References 



Quantifying spatially heterogeneous dynamics in computer simulations of glass-forming liquidsll 



M. D. Ediger, C. A. Angell, and S. R. Nagel 1996 J. Phys. Chem. 100 13200. 

W. Gotze and L. Sjogren 1996 Chem. Phys. 212 47; W. Gotze and L. Sjogren 1995 Transport 

Theory Stat. Phys. 24 801. 
H.Z. Cummins, these proceedings. 

See the special issue of Transport Theory Stat. Phys. 24 devoted to experimental and 

computational tests of MCT. 
See, e.g., A. van Blaaderen and P. Wiltzius 1995 Science 270 1177, and R. Leheny, et al 1996 J. 

Chem. Phys. 105 7783. 
K. Schmidt-Rohr and H.W. Spiess 1991 Phys. Rev. Lett. 66 3020. 

R. Bohmer, R.V. Chamberlain, G. Diezemann, B. Geil, A. Heuer, G. Hinze, S.C. Kuebler, R. 

Richert, B. Schiener, H. Sillescu, H.W. Spiess, U. Tracht and M. Wilhelm 1998 J. Non-Cryst. 

Solids 235-237 1, and the many references therein. 
See, e.g., M.T. Cicerone and M.D. Ediger 1995 J. Chem. Phys. 103 5684. 
For a review, see H. Sillescu J. Non-Cryst. Solids in press. 

W. Kob, C. Donati, S.J. Plimpton, P.H. Poole and S.C. Glotzcr 1997 Phys. Rev. Lett. 79 2827. 
C. Donati, J.F. Douglas, W. Kob, S.J. Plimpton, P.H. Poole and S.C. Glotzer 1998 Phys. Rev. 
Lett. 80 2338. 

C. Donati, S.C. Glotzcr, P.H. Poole, W. Kob, and S.J. Plimpton, |cond-mat/9810060| . 
S.C. Glotzcr, C. Donati and P.H. Poole 1998 in Computer Simulation Studies in Condensed- 
Matter Physics XL, edited by D.P. Landau, ct al. (Springer- Verlag, Berlin) (in press). 
P.H. Poole, C. Donati and S.C. Glotzer Physica A in press. 



C. Donati, P.H. Poole and S.C. Glotzer, cond-mat/9811145 



P. Allegrini, J.F. Douglas and S.C. Glotzer, in preparation; P. Allegrini and S.C. Glotzer, in 
preparation. 

C. Bennemann, C. Donati, J. Baschnagel, and S.C. Glotzer, preprint. 

J. P. Hansen and I.R. McDonald 1986 Theory of Simple Liquids (Academic Press, London). 
H.E. Stanley 1971 Introduction to Phase Transitions and Critical Phenomena (Oxford University 
Press, New York). 

A related correlation function has been calculated for the Ising spin glass in S.C. Glotzer, N. Jan, 

T. Lookman, A.B. Maclsaac, and P.H. Poole 1998 Phys. Rev. E 57 7350. 
The Lennard- Jones interaction parameters e a $ and o~ a p are given by: caa = 1-0, cab = 1-5, 
ess = 0.5, oaa = 1-0, <jab = 0.8, obb — 0.88. Lengths are defined in units of oaa, 
temperature T in units of (AA/ks, and time t in units of a\ A m/ eAA- Both types of particles 
are taken to have the same mass. 
W. Kob and H. C. Andersen 1994 Phys. Rev. Lett. 73 1376; W. Kob and H. C. Andersen 1995 

Phys. Rev. E 51 4626; 1995 Phys. Rev. E 52, 4134. 
Correlations of particle motions on short time scales has been studied in 2-d simulations by Y. 

Hiwatari and T. Muranaka 1998 J. N on- crystalline Solids 235-237, 19, and references therein. 
Indirect evidence for a growing length from studies of glassforming liquids or polymers in confined 
geometries is reported in, e.g. M. Arndt, et al., Phys. Rev. Lett. 79, 2077 (1997); B. Jerome 
and J. Commandeur, Nature 387, 589 (1997); and P. Ray and K. Binder, Europhysics Lett. 27, 
53 (1994). A "hydrodynamic" correlation length was calculated for a LJ liquid by R. Mountain 
in Supercooled Liquids: Advances and Novel Applications, ACS Symposium Series 676, 122 
(1987). 

[25] P. Allegrini, C. Donati, J.F. Douglas and S.C. Glotzer, in preparation. 

[26] D. Stauffer 1995 Introduction to Percolation Theory (Taylor & Francis Ltd, London). 

[27] Other recent computational investigations of spatially heterogeneous dynamics include Y. 

Yamamoto and A. Onuki 1998 Phys. Rev. E 58 3515; D. Perera and P. Harrowell 1998 J. 

Non-cryst. Sol. 235-237 314; B. Doliwa and A. Heuer 1998 Phys. Rev. Lett. 80 4915; Ref. §3l 



